home *** CD-ROM | disk | FTP | other *** search
/ Acorn User: China / Acorn User China CD-ROM (UK) (Disc B) / Acorn User China CD-ROM (UK) (Disc B).bin / STUTTGART / LANG / GOFER / scripts_InterPol < prev    next >
Encoding:
Text File  |  1993-11-24  |  1.4 KB  |  42 lines
  1. -- Interpolation of integer functions by polynomials
  2. -- with integer coefficients.
  3.  
  4. -- coeffs of polynomial of degree n agreeing with f at 0,1,..n
  5. coeffs :: Int -> (Int -> Int) -> [Int]
  6. coeffs n f = [ coeff_to n k f | k <- [0..n] ]
  7.  
  8. -- for example Go> coeffs 6 (\x->(x^2-1)^3)
  9.  
  10. -- Convolve with stirling numbers to get k-th coefficient of function f,
  11. -- assuming a polynomial of degree n.
  12. coeff_to :: Int -> Int -> (Int -> Int) -> Int
  13. coeff_to n k f = alt_sum [s*b | (s,b) <- zip ss bs]
  14.         where ss             = [ stir1 (k+i) k | i <- [0..n-k] ]
  15.               bs             = drop k (take (n+1) (diffs f))
  16.               alt_sum []     = 0
  17.               alt_sum (x:xs) = x - alt_sum xs
  18.  
  19. -- stir1   Stirling numbers of the 1-st kind
  20. stir1 :: Int -> Int -> Int
  21. stir1 n k | k < 0     = 0
  22.           | n < k     = 0
  23.           | n == k    = 1
  24.           | otherwise = (n-1)*(stir1 (n-1) k) + stir1 (n-1) (k-1)
  25.  
  26. -- Get the Newton series of iterated differences of a function f
  27. diffs :: (Int -> Int) -> [Int]
  28. diffs f = [ (u 0)/v | (u,v,_) <- iterate step (f,1,1) ]
  29.             where step (h,d,n) = (delta h,n*d,n+1)
  30.                   delta h x    = (h (x+1)) - (h x)
  31.                   iterate g y  = y:iterate g (g y)
  32.  
  33. zip :: [a] -> [b] -> [(a,b)]
  34. zip as []         = []
  35. zip [] bs         = []
  36. zip (a:as) (b:bs) = (a,b):zip as bs
  37.  
  38. drop :: Int -> [a] -> [a]
  39. drop 0 xs         = xs
  40. drop _ []         = []
  41. drop (n+1) (x:xs) = drop n xs
  42.